Nonlinear optical response in gapped graphene 
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We present a formulation for the nonlinear optical response in gapped graphene, where the low-energy single- 
particle spectrum is modeled by massive Dirac theory. As a representative example of the formulation presented 
here, we obtain closed form formula for the third harmonic generation (THG) in gapped graphene. It turns out 
that the covariant form of the low-energy theory gives rise to a peculiar logarithmic singularities in the nonlinear 
optical spectra. The universal functional dependence of the response function on dimension-less quantities 
indicates that the optical nonlinearity can be largely enhanced by tuning the gap to smaller values. 
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INTRODUCTION 

Graphene is the first example of the realization of a truly 
two-dimensional (2D) crystal made of carbon atoms HI]. This 
intriguing material in addition to promise for novel applica- 
tions, from a fundamental point of view, provides the con- 
densed matter community with a low-energy laboratory of 
Dirac electrons on the table top 101 ■ Since then there has been 
remarkable success in pushing the idea of employing elec- 
tronic, mechanical and various other properties of graphene in 
technological applications IslSl- The chiral nature of Dirac 
charge careers in graphene makes them non-stoppable which 
mathematically manifests as the absence of back-scattering. 
This means that with pristine graphene which contains mass- 
less Dirac fermions, there can be no "off state in electronic 
applications. Therefore, to enable the use of graphene in elec- 
tronic devices, one needs to open up a gap in the single particle 
spectrum jst] usually by statically reducing the symmetry via 
extrinsic effects 101, or by coupling to another field to gener- 
ate dynamical masses 101 ■ Massive Dirac fermions possess a 
single-particle gap which causes the graphene to behave like 
a truly 2D semiconductor in some respects. Nevertheless the 
nature of such "relativistic" massive theory is drastically dif- 
ferent from the usual parabolic bands in a semi-conductor 

Before the synthesis of graphene, the problem of two- 
dimensions has been usually approached from the third di- 
mension by e.g. a geometrical confinement in hetero- 
structures [SD, or by appropriately chosen B-field which ef- 
fectively confines the dynamics of carriers into two spatial di- 
mensions In this respect, truly 2D gapped graphene pro- 
vides a novel platform for non-linear optical applications. The 
optical properties of 2D graphene have been extensively inves- 
tigated at linear level lid, ll 111 . However, at nonlinear level, 
there has been a limited number of studies: A classical theory 
for the electromagnetic (EM) response of pristine graphene 
was developed by Mikhailov and Zeigler 11211 who attributed 
strong nonlinear EM response to the massless nature of Dirac 
electrons in graphene. Nonlinear current response of massless 



'Electronic address: akbar.jafari@gmail.com 



Dirac fermions was studied by Wright and coworkers, who 
used direct expansion of the wave function in terms of mul- 
tiples of frequency of the applied electric field 11311 . They 
found a high triple-frequency current response in massless 
graphene. On the experimental side, Hendry and collabora- 
tors used four-wave mixing to study the nonlinear optical re- 
sponse in graphene flaks, where they found a remarkable large 
third-order optical response lfl4ll . Also broadband optical non- 
linearity was observed by Wang and coworkers in graphene 
dispersions 1 15]. 

Therefore it is timely to consider the problem of third order 
optical response in the 2D lattice of graphene by generaliz- 
ing it to the massive case. The problem of third order op- 
tical response for l-i-l dimensional massive Dirac fermions 
has been previously considered by Wu flit. So the present 
work can also be considered as a natural generalization of the 
Wu's work to 2-1-1 dimensions. In one spatial dimension, or 
for quasi one-dimensional systems such as organic materials 
or carbon nano-tubes closed form expressions for the optical 
non-linear responses can be found llTlllSll . We employ a for- 
mulation we have developed earlier for investigation of the 
nonlinear optical properties in the matrix form 12011 . We find 
that the covariant form of the dispersion relation for massive 
Dirac fermions in 2-1-1 also allows for a simple closed form 
expression for the third order response in general. We de- 
velop the general theory theory of nonlinear optical response 
for 2-1-1 Dirac fermions with arbitrary gap parameter in terms 
of the Feynman diagrams. Then as a representative exam- 
ple, we evaluate the line-shape of the third harmonic genera- 
tion (THG). We obtain a universal logarithmic functional form 
which depends only on the combination ly/m, where i/ is the 
photon energy, and m is the gap parameter. 



MODEL AND METHOD 

The single-particle energy band structure of Graphene con- 
sists of two Dirac cones which are connected to each other 
by time-reversal symmetry. For optical applications the mo- 
mentum of light compared can be safely ignored and hence 
all optical processes of interest will take place around a single 
cone. Therefore it is sufficient to consider only one valley cor- 
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responding to which a spinor -ipt = (at, bt) denotes creation 
of pz electrons of momentum p in sub-lattices A, B, respec- 
tively. Then the effective low-energy Hamiltonian around the 
valley under consideration can be written as, 



H = ^ -ipivpip.a + ■mvF(Jz)i^p 



(1) 



where p ~ Px + ipy = P ^^^'^ is the complex representa- 
tion of a 2D vector. For simplicity we work in units with 
vp = e = = ^- The physical constants can be restored at 
the end of calculations if required lll9ll . The above Hamilto- 
nian can be brought into diagonal form if the following unitary 
transformation is applied 




2£p \ — m —y/Sff+ m 



to induce the change of basis, 



hp 



(2) 



(3) 



to a new basis of conduction (Cfj) and valence (vp) bands, (pp = 
{cp, vp), with correspondingly positive and negative energies 
Cp = ±sp. The diagonal form of the Hamiltonian in the new 
basis, </)p- reads. 



(4) 



To proceed further, we need the explicit expression for the 
current operator of 2h-1 dimensional Dirac electrons in the 
new basis. Therefore we derive in the following the current 
operator as follows: The current operator in terms of original 
real space spinor ipp can be written as. 



J 



E4 



d_ 

dp 



[p.a + m(Tz) 



The transformation ipp = Vp4)p on the spinors, gives the fol- 
lowing form. 



■h = ^ 4>pVpaiVp4>p, i 



x,y 



(6) 



where the transformed Pauli matrices are given by, 

. -mcos(/3j? pcosifp 
VJ. (7,r Vff = + sm tpp ay H (7) 



'p^x "p 



y^p^yVp 



ep 

-msaupp 



p sm ipp 
- cos ipp Oy H (Tz (8) 



of the vector character of the three Pauli matrices under SO(3) 
rotations. 

MULTI-CURRENT CORRELATIONS 



As detailed in Ref. |20ll . we are interested in calculation 
of multi-current loops to which photon propagators are at- 
tached. As a demonstration of this method, in the following 
we first calculate the two-current response which is connected 
to optical conductivity. After reproducing well-known results, 
along the same lines we proceed to calculate four-current cor- 
relations within our matrix diagrammatic formulation which 
is suitable for situations such as massive Dirac fermions in 
gapped graphene. 



zi = lujn + ifi 




23 = IWn + IVa 



FIG. 1: A typical Feynman diagram corresponding to a four- 
operator correlation function (ABCQ) contributing in the nonlinear 
optical response. Operators A, B, C, Q could be any matter operator 
that couples to some power of the gauge field A of the incident light 
denoted by wavy lines (photons). The frequency Va is the sum of all 
incoming frequencies vi + 1^2 + va- 



in general the fully retarded expectation value of nested 
commutators of various operators, e.g. A, B, C, Q schemati- 
cally depicted in Fig.[T]can be Lehman represented as 1 23, 21 1, 



AoaBabCbcQ 



cO 



abc V 



(9) 

where the sum of frequencies is — vi + V2 + 1^3, and 
J2-P permutes {A, vi), {B, (C, 1^3), (Q, -i^^) around the 
current loop, and Bab, etc. stand for the matrix element 
{a\B\b). Note that in the above formula, the substitution 
— > = J/ + iO for all frequencies is understood. 



Linear response: Optical conductivity 



Note that the second equation can be obtained from the first 
one by simply replacing ip (p—TT/2, which is a consequence 



Within this formulation, the two-current response is given 
by ii. 
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1 1 



-Gx + sm ipp ay 



-mcosff; 



■ sm 03^ (Jr, 



(10) 



where iv is the external (photon) frequency. The correspond- 
ing diagram is simpler than the one in Fig. [T] where only 
two photons are attached, and only two frequencies Za = 
iuJn+iav, a = 0, 1 are present. In the above equation, quanti- 
ties in brackets are inter-band terms of the current operator and 
contain all necessary information about the matrix elements of 
the current operator between the valence and conduction band. 
The matrix multiplications required in the above expression 
are easily performed to simplify the result as. 



(J.J.)=2!Tri-l^^ 

(11) 

Note that the sum over permutations of two Jx operators 
around the loop produces a 2! factor in the right hand side ll20ll . 
as the above formula is symmetric with respect to fermion 
propagation frequencies zq, zi. Standard Matsubara summa- 
tion, the two-current correlation function simplifies to. 
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(12) 



where the temperature is assumed to be zero so that the Fermi 
function is 1 only for negative energies, and is zero otherwise. 
After analytic continuation, iv v + iij, the spectral function 
(imaginary part) corresponding to the absorption - i.e. the first 
term in the bracket - can be easily obtained. To compare with 
the existing results in the literature, we use the conductivity 
formula (Jxxii^) = ^"^{JxJx) I ij-v) and restore the fundamen- 
tal constants to obtain 



K&Gxx{v) 



e 

8^ 



1 



4m^ 



2 



1 



(13) 



which after noting the fact that our calculation has been done 
for a single-valley, agrees with e.g. eq. (10) of Ref. 12211 . Note 



that here = j{^K) is the conductance of ideal graphene. 

In a similar way, one can calculate correlation functions like 
{JxJy)- In simplifying such elements one should note that 
odd powers of sin tp or cos ip average to zero upon ip integra- 
tion. Moreover, summation over permutations V picks up the 
most symmetric part. In the case of two-current correlation 
the xy component will become zero, unless a magnetic field 
is apphed 1123 1. 



Four-current correlations 

Now that we checked our formulation with well-known re- 
sults, let us calculate higher order responses. But before pro- 
ceeding to question of higher order responses, we note that 



in a general material the higher order correlation between the 
current operators, i.e. {JaJbJcJd) which essentially contains 
four velocity vertices along a, fo, c, d G {x, t/} directions of 
space, is not the only important term when one is interested in 
third order response. In general terms of the type JaTbcJd, etc 
also may appear where the stress tensor t^c too contributes to 
the third order response 12 ill . Similar to current vertices given 
by 



even-parity vertices of the form 

dpcdpd 



(14) 



(15) 



may also enter the theory of nonlinear optical response for 
materials with arbitrary band structure defined by ep. In the 
case of ideal graphene (m = 0), it is interesting to note that, 
due to the linear form of energy as a function of momentum 
p, only velocity vertices will appear in the theory as we cal- 
culated them in (|6|. All higher order derivatives related to the 
stress tensor, etc. will be identically zero for massless Dirac 
fermions. In the case of massive Dirac fermions where the 
asymptotic behavior of the energy dispersion for photon en- 
ergy scales much higher than the gap parameter m, becomes 
linear, the terms arising from stress tensor will have negligi- 
ble contributions, for energy scales much beyond those corre- 
sponding to band edge excitations. 

With this point in mind, in this paper we focus on four- 
current correlations of the form {JxJxJxJx) ■ From the right- 
most current operator, only the inter-band terms contribute 
which lead to an excitation into the conduction band. Then 
the second operator from right will have to create intra-band 
excitation. However, for the third and fourth operator, there 
are two different possibilities depicted in Fig. |2l Either the 
operator number 3 creates an intra-band excitation as in part 
(i) of Fig. ID In this case the operator number 4 must create 
an inter-band excitation. Such sequence of operators can be 
represented by — h H — , where — denotes /nfer-band vertex, 
and + stands for an /«fra-band vertex. With this notation, the 
second possibility denoted in part (ii) of Fig. |2]can be sum- 
marized as H 1 — , i.e. the third operator returns the sys- 
tem to valence band, and operator number 4 has to create an 

intra-band excitation. Note that a sequence like is 

not a genuine four-operator correlation, but rather a product 
of two-operator correlations and can be extracted from lin- 
ear responses as well. All other possibilities can be generated 
by the sum over permutation X^p- Therefore the above se- 
quences are two general families of excitation patterns which 



4 




matrix. Adding the above terms gives. 



p cos ( 
p 



el 
p 



[zqZi + Z0Z2 - Z1Z2 - el] 



X 2el X 



(18) 



Now we can evaluate the Matsubara sums, and at the end the 
resuh will be symmetrized with respect to exchange of the 
frequencies of the attached photon vertices. Using the angular 
averages, 



FIG. 2: (Color online) Schematic sketch of two types of operator se- 
quences contributing to the four-operator correlation function. Type 
(i) and (ii) processes are denoted by green and blue, respectively. 
Indices 1 to 4 label the operators from right to left. 



(COS^ (fi) 



6 



(cos^ If sin^ cp) 



the angular integration can be simplified to, 



1 



(19) 



by the cyclic property of the trace involved in the loop diagram 
of Fig. [1] include all possibilities. These are summarized in 
Fig.E] 

Now let us calculate the contribution of two sequences of 
excitations (i) and (ii). First note that the only terms surviving 
the trace operation are those proportional to unit 2x2 ma- 
trix, and that odd powers of sin ip or cos (p do not survive the 
angular integration arising from J^p- Therefore, the non-zero 
terms of the first class of terms can be simplified after some 
algebra as. 



1 



3 I 

II ^2 _ p2 

.a=0 " P. 



(16) 



2 2 
p cos I 

el~ 

^p 



m} cos^ pf 2 
p V sin kpp 



[{ziZ2 + el){zoZ3 - el) + el{zi + Z2)(^o - 23)] , 

where frequencies Zq, a = 0, 1, 2, 3 are defined in Fig. [T] 
Similarly for the class (ii) terms we obtain. 



<Pp 
J (27r)2 



3 1 

II 72 _ p2 

.a=0 " ^P. 



(17) 



p2 COS^ Iy9p- 

el ■ 

^p 



m? cos^ ipp 2 

-p h sill ipj! 



[£%Z2 - Z'i){zQ - zi) - {el - Z2Z'i){el - zazi)] , 



where an overall factor of 2 arises from trace over the unit 



^ J 27r 



V 

2 



el 
p. 
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et 
p 



n 

a;=0 ' 

[zoZl + Z0Z2- Z1Z2 



4] 



(20) 



Up to this point, the formulation is quite general. At this point, 
let us specialize to the case of THG, where wi ~ iv2 = iv^ = 
iv and ivc ~ Ziv. In this case, the Feynman diagram shown in 
Fig.[T]must be labeled with z^ ~ iojn + aii' for a = 0, 1, 2, 3. 
Permuting the external vertices amounts to moving the out- 
going 3iiy photon frequency around the loop. Such permuta- 
tion can be accounted for by a cyclic permutation of the set 
{^0, zi, Z2, Z3}. Performing the J^-p will only affect the last 
term in Eq. (|20|) as. 



^02^2 



Z1Z2 ~ el] =2 [zqZ2 + Z1Z3- 2el] 



where the THG identification, z„ = ia;„ + aii/ with a ~ 
0, 1, 2, 3 is understood. Let us emphasize that the above result 
for the sum over permutation holds only for THG. Therefore 
the THG susceptibility will become. 



THG 



f Pdp XT- 1 



m 

6^ + 1 

et 

p 



el 

Q = " P. 



[zoZ2 + Z1Z3- 2el] 



(21) 



The 1 //3 performed with standard contour inte- 

gration techniques producing Fermi functions /{ztep) corre- 
sponding to poles at conduction and valence bands, respec- 
tively. Assuming the temperature to be zero, and for undoped 
massive Dirac spectrum, only contributions from poles at the 
valence band will be left for which the residues can be calcu- 



5 



THG in Gapped Graphene 




FIG. 3: (Color online) THG line-shape for massive Dirac fermions in 
gapped graphene. The gap parameter m is taken to be unit of energy 
in this figure. Fundamental constants are assumed to be e = /t = 1. 



lated to give. 



.2 

4e|- 



(22) 



2 ep{v'] 



94) 



Note that the analytic continuation iv ^ = v 
been performed in the above equation ll20ll . 



iO has 



RESULT AND DISCUSSIONS 



The radial momentum integration in Eq. (l22l i can be per- 



formed with the change of variable to 
gives the following closed form result 



7,^ which 



X 



THG 



(^+) 



1 
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+ {Arr? - 9i/+)(8m2 + 3i^+) In 
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2m 


+ 'iv+ 


2m 


— 3i^+ 



(23) 



First of all, note that the above expression when the funda- 
mental constants are restored will be proportional to jY? . 
Therefore the above expression shows the results in the nat- 
ural units. A remarkable feature of the above expression is 
that, the quantity my^^'^ is only a function of vjm. This 
functional form is universal characteristic of 2h-1 dimensional 
Dirac fermions. Hence in Fig.|3] we have plotted the real part 
and the magnitude of THG line-shape incorporating this ob- 
servation. This peculiar scaling form indicates that in the limit 



where m — > 0, and the system approaches the gapless limit, 
since mx^^*^ in the natural units employed here will have to 
remain on the scale of unity, the x^^^ itself will grow in- 
versely proportional to the gap parameter m. This sheds a 
new light into why pristine (gapless) graphene is expected to 
display large optical nonlinearity lll4ll . 



Apart from the zero frequency divergence, which is a gen- 
eral feature of optical response functions, an interesting point 
in the above expression which can be noticed is that, it is pre- 
cisely the covariant (relativistic) form of the dispersion rela- 



tion £ 



+ m^, which after changing integration variables 



from momentum p, to energy e leads to logarithmic depen- 
dence at finite frequencies f^/m ~ 2/1, where £ is an in- 
teger £ — 1,2,3. This would not be the case for a typical 
parabolic dispersion relation in ordinary semiconductors. The 
main THG peak corresponding to £ = 1 and the second one 
corresponding to ^ = 2 are clearly seen in the THG spectrum 
depicted in Fig. [3] The feature corresponding to ^ = 3 is 
the weakest feature. Note that as we calculated in Eq. ( fT3] ), 
such logarithmic dependence does not show up in the linear 
optical response. Therefore it can be considered as a valuable 
information which is contained only in higher order optical 
responses. Moreover, this logarithmic dependence is not re- 
stricted to the THG spectrum, but it also appears in all other 
higher order responses. To see this, we note in Eq. ( l20b that 
in the special case of THG, only the functional form of 
as a function of loop frequencies Zq, will be affected. But 
the dimensional form of the summand will generally retain a 
similar form. 

The calculation of THG for massive Dirac fermions in 1h-1 
dimension by Wu 1161 indicates a characteristic inverse square 
root divergence at the main THG peak v^^i and its harmon- 
ics. In the case of small clusters (of correlated 2D nature), 
a numerical investigation of THG spectrum suggests a line- 
shape composed of a superposition of resonances correspond- 
ing to simple poles 12411 . For extended 2D system with broken 
symmetry, it was shown that the nesting underlying spin den- 
sity wave instability can produce a peculiar z~'^/'^ lii(z) form 
of singularity in the non-linear optical spectra |20l]. With the 
above examples in mind, the gapped graphene as a first exam- 
ple of truly two dimensional system provides us with a unique 
example of a 2D system where optical nonlinearity mani- 
fests in a universal logarithmic function of the dimension-less 
quantity v/m. 

It is interesting to check the limit of gapless graphene by 
taking the limit to ^ in Eq. ( |23] ). which gives, 



mx 



THG 



54i TO 
^ 192^ 



(24) 



The coefficient of the above expression is universal and can 
be regarded as a hallmark of the ideal graphene at the THG 
spectroscopy. 
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